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In this paper we use different techniques from the fractional and pseudo-operators calculus to 
solve partial differential equations involving operators with non integer exponents. We apply the 
method to equations resembling generalizations of the heat equations and discuss the possibility of 
extending the procedure to the relativistic Schrodinger and Dirac equations. 



I. INTRODUCTION 

Evolution equations involving fractional derivatives play an important role in the theory of transport in heteroge- 
neous media [1] and different techniques have been developed [2], or are under study [3], to deal, in an unambiguous 
way, with the underlying formalism. 

This paper is a contribution to the theory of operators involving noninteger exponents and we will show how many 
of the techniques employed within such a context can be used to get further insight into the formalism of differential 
equations, like the relativistic Schrodinger equation or the pscudohcat equation. Here we will be mainly concerned with 
the mathematical aspects and postpone the application of the underlying formalism to a forthcoming investigation. 

Different mathematical tools to deal with fractional differential equations are available. Many of them were not 
considered by physicists during the second decade of the last century, at the time of the proposal of the Dirac equation. 
Presently some progress has been made toward the understanding of fractional operators and their use is not limited 
to pure mathematics anymore. However, it seems that a good convergence can be realized by at least three different 
points of view, namely the formalism of pseudodifferential operators, the fractional calculus and the theory of integral 
operators, as will be shown in this paper. 

We will discuss evolution equations containing the square roots of differential operators and we will adopt different 
methods of solution, involving the techniques from fractional calculus, from the theory of pseudo-operators [4], and 
from that of the generalized integral transform [5]. In the following we will see that the use of a judicious combination 
of operational techniques and of integral transform methods is one of the main ingredients to get a solution of partial 
differential equation involving fractional derivatives. 

Evolution equations expressed in terms of the square root of an operator have a quite old story in physics and an 
example is provided by the debate developed around the meaning of the relativistic Schrodinger equation, which paved 
the way to the formulation of the Dirac equation for the electron [6] . As a first example of this class of equations, we 
consider, without mentioning the specific physical problems it describes, the following equation (the x and r variables 
are dimensionless): 

drF{x, r) = - dl'^ F{x, t) , F{x, 0) = f{x) . (1.1) 

The relevant solution can formally be written as 

F(x,r) = exp(-Tay^) /(x), (1.2) 

that, unless we specify the action of the exponential operator containing the square root of a derivative on the initial 
function f{x), is just a way of restating Eq. (1.1). The Doetsch transform [7] (see Appendix A) is the appropriate 
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tool to get a solution for our problem, that can be formally written as 
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and is valid if the integral on the right-hand side converges. As we will see in the following, the fact that the solution 
of Eq. (1.1) is given by the form of an integral representation is a common feature of problems involving fractional 
derivatives [8]. 

We consider now the equation 
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to which we will refer as the pseudofractional (or simply pseudo) heat equation. By applying the evolution operator 
formalism, the solution of this equation can be written as follows 



Fix,T) = e-^V^.fix), 
and, according to the Doetsch transform, we get 
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The crucial step to achieve the complete solution consists in specifying the action of the exponential involving the 
second-order derivative. This can be done by using the Gaiiss transform [9] 



e"^^/(x) = r d^expj-^^--^^ /(O, 



that allows us to write 
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As a particular case, we note that if f{x) = e we can use the so called Glaisher identity [10] 
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that, once inserted in Eq. (1.6), provides the result 
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As clearly shown in Fig. 1, compared to the ordinary heat diffusion equation, the pseudoheat equation yields an 

evolution characterized by a less significant spreading and a larger reduction of the peak. 

More conventional methods, based, for example, on Fourier transform techniques can also be used. If we assume 
that the function f{x) admits a Fourier transform, namely 



m - 

the solution of our problem can be obtained as^ 
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^ In the second line we have used the identity f{dx) e" ^ = /(a) e'^ 



FIG. 1; The evolution of a packet initially Gaussian (blue) according to an ordinary (green) and a fractional (purple) heat 
equation. The case considered corresponds to r = 1. 

Solutions (1.8) and (1.12) are both expressed in terms of an integral transform and, in principle, there is no reason to 
prefer one over the other except for the nature of the initial function and the convergence of the integral involved in 
the expressions. 

The use of the Fourier transform method is a useful alternative to the method we are discussing when the pseudo- 
operator is not a function of the derivative operator only. To clarify this point we note that the evolution equation 

dr F{x, t) = P(a,) F[x, t) , F{x, 0) = f{x) , (1.13) 

with P(.) an operator or an analytical function, can be solved as 

F{x,t) ^ dfce^-f'(*'^)/(fc)e*'=^, (1.14) 

V27r J^oo 

which holds if the integral converges. On the contrary, if P(.) is a function of the derivative operator and the 
coordinate, this method, even though still applicable (at the price of using a more complicated family of "symbols"), 
is more cumbersome and an alternative procedure might be more convenient. This is the case of the following equation 

drF{x,T) = -^x - cd,F{x,T) , F(x,0) = fix), (1.15) 

where c is a constant. Here, the use of Doetsch transform and Weyl disentanglement formula^ allow us to write the 
relevant solution in a quite straightforward way as follows 

F(a:,r) = -i- H dt cxp \ - ^- tA - ^^^^ f{x) (1.16) 



= f ^« «p { - ( + ' + ^) } «- - /(- + , 

while the use of pseudodifferential operators would require a more cumbersome analysis. 

According to the discussion developed so far it is evident that the combined use of integral transform and operational 
methods offers the natural environment to deal with fractional and/or pseudo-operators. Within this context the 
Laplace transform method offers a wealth of possibilities. For example, by assuming that the identity [7] 
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can be extended to operators, we find [11] 



1 r°° 

== = / dfJoWe-*^^. (1.18) 
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We can therefore express the solution of a problem of the type 

+ lf{x) = g{x), (1.19) 

with f{x) an unknown function, as 

fix) = / dtJo{t)g{x - t) . (1.20) 
Jo 

The use of the above transform allows to cast the equation 

drF{x,T) = ^dl + lF(x,r) , (1.21) 

in the form 

/•oo 

F(x,t) = \ dtJo{t)drF{x-t,T). (1.22) 
Jo 

Unfortunately, this equation is of limited usefulness since the convolution integral on right-hand side does not produce 
well behaved functions. 

II. RELATIVISTIC SCHRODINGER EQUATION 

As a variation on the pseudo-fractional heat equation, we consider the relativistic (1+1 dimensional) Schrodinger 
equation for a particle of mass m, that writes 

idr^{v,T) = ^1 - d^^{r],T) , *(7?,0) = ^iv), (2.1) 

where time and position variables are normalized to the Compton wavelength (^c = h/mc): 

ct X 

It describes the quantum evolution of a relativistic free particle and does not include negative energy contributions. 
The use of the Schrodinger equation in the form given by Eq. (2.1) has been criticized on the basis of arguments 
concerning its nonlocality, as a consequence of the nonlocal character of the fractional operators. The problem is of 
genuine mathematical nature and does not arise from the potential as in the case of the nonlocalities associated with 
the fermionic nature of nucleons in nuclear many-body problems. The initial negative reaction to Eq. (2.1) stems, 
probably, from the intrinsic mathematical difficulties rather than due to a real physical reason. In Ref. [12] it is 
clearly stated that Eq. (2.1) has been ruled out as a tool to develop the relativistic quantum theory for reason of 
simplicity (i.e., to deal with simpler expressions, avoiding square roots of operators) but not on the basis of a real 
"physical cogency." Even though the questions about the intrinsic nonlocal nature of Eq. (2.1) are subtle and have 
opened a long-standing debate, the relevant physical consequences do not display any disturbing feature. This point 
has been reconsidered in various papers (see Ref. [13]) where its physical legitimacy has been discussed in depth. 

From the mathematical point Eq. (2.1) can be transformed into a pscnidoheat equation by a Wick rotation, and, 
thus, its solution can be written by following the procedure outlined in Sec. I. However, the obtained result is 
ambiguous and hampered by the relevant physical interpretation and the fact that the convergence of the integral 
representation (1.8) with — is not ensured. For this reason we will treat the problem following a different 

approach. 

The evolution operator associated to Eq. (2.1) can be written as follows 
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The fractional operator appearing in the series can be handled by means of a trick, often exploited in the theory of 
the fractional derivatives (see also the forthcoming section), namely 



^'"-m, (2.3) 



I 1 - 1 - r°° e"* 

where the last equality has been obtained using the Laplace transform identity 
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In the case (^(r/) = e , from the Glaisher identity (1.9) we obtain 

yrr^e-''^ = ^5 r d. ^ ' exp(- + -^n , 
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and taking into account the following identity involving the two- variable Hermite polynomials [11, 14] 



(2.5) 
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we can cast the free particle solution of the relativistic Schrodinger equation in the form 

*(?7, r) = A{r,, r) e" + i Bin, r), (2.7) 
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The comparison of the previous solution with its non-rclativistic counterpart shows significant differences, which 
becomes more evident with increasing time (see Fig. 2). We will discuss the relevant physical meaning in the 
forthcoming sections. 

Even though the above relation provides us with a well behaved solution, extendable to non-Gaussian packets (see 
the concluding section), its unappealing features stems from the fact that it is essentially a Taylor series expansion in 
time. To overcome this problem different strategies (including the already quoted Fourier transform method) will be 
developed in the forthcoming sections. 

Another commonly adopted method in fractional calculus, which allows the elimination of noninteger exponents, is 
based on the elementary Laplace transform identity already exploited to derive Eq. (2.5). By using this identity, we 
can recast Eq. (2.1) in the form 

idr^{v,T) = -d^^^r,,T), (2.9) 

where 

*(??, r) = exp - i , I *(??, r) 

$(77,r) = ^ / ds^e"^^ *(r/,T) (2.10) 
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FIG. 2: Time evolution of an initially Gaussian packet (blue) undergoing diffusion ruled by a relativistic Schrodinger equation: 
r = 0.5 (green) and r = 1 (purple). 
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FIG. 3: The function $ of Eq. (2.10) (green) calculated for an initial wave function *(C, 0) = e" ^ (blue). 

This is a Schrodinger equation in which time and spatial derivative operators do not act on the same wave function. 
The action of the second-order spatial derivative is indeed mediated by another operator, inducing a kind of diffusion, 
whose effect is that of delocalizing the initial wave function, as is illustrated in Fig. 3, where we have reported how 
the integral transform in the previous equation affects a given function. An alternative way of writing Eq. (2.9) is 
the following 

z5,*(7?,r) = -52^^'(ry,T), (2.11) 

with 

D = ^ d.^e^^., (2.12) 

that allows the series solution 

*(^'^) = E ^2.13) 

n=0 
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where 



= £)§„_i(r?) = J- ds ^ d^ exp |- a| . (2.14) 



Let us now consider the inclusion of a potential, so that Eq. (2.1) becomes 

idr^{V,r) = yi - 52 + V{v)] ^V,r). 
In this case, by introducing the transformation (a kind of interaction picture) 

A(ry,T) = e-'^V^^^{r],T), 

the problem to solve becomes 

idrA{r],T) = V{r],T)A{r],T) 

where 



(2.15) 

(2.16) 
(2.17) 



The previously outlined procedures can still be applied, but the solution of problem (2.17) should be expressed as it 
follows: 



*(ry, r) = e* " | exp i £ dr' V{x, r')^ | *(r?, 0) , 



(2.18) 



where {}+ denotes the Dyson time-ordering operation. In this case the ordering procedure becomes mandatory 
because the Hamiltonian contains the explicitly time-dependent term V{r],T), which does not commute with itself at 
different times. 



III. THE HEISENBERG PICTURE AND THE FRACTIONAL OPERATORS 

In this section we will show how the use of methods involving fractional derivatives can be useful to study physical 
states whose evolution is ruled by relativistic Hamiltonians. We will show that these methods find a natural application 
in the treatment of the Heisenberg equations of motion of the physical observables and in the evaluation of their average 
values. 

The relativistic Hamiltonian for a free particle is (we limit ourselves to the one-dimensional motion) 

H = c + p^, (3.1) 

and the Heisenberg equations of motion for the position and momentum operators can be written as^ 

x,H 



d . _ 1 
dt ih 



^ ^ ^P = 0, (3.2) 

\/w? (? + fP' dt 



and the relevant solution reads 



m = m^^ "fl t. (3.3) 

As also noted in Ref. [13], the operator dxjdt is the velocity operator, whose definition has been introduced in quite 
a natural way within the present context, without the necessity of defining a new position operator as usually done 
in the case of the Dirac theory and as discussed in the forthcoming section. 



3 Note that [£,/(j5)] = ihdpf{p). 
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The dynamical behavior of a wave packet undergoing an evolution ruled by the Hamiltonian (3.1) can be inferred 
from the evaluation of the particle position and momentum at a given time t. By assuming a packet initially Gaussian 
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we find that it spreads in time according to the relation 
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The evaluation of the contribution depending on the momentum operator can be performed using identities (2.4) and 
(1.9), and yields 



(3.6) 
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The explicit evaluation of the previous integral allows to write the width of the packet in the form 



where 
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Equation (3.7) reduces to its nonrelativistic version for R{a) = 1. As shown in Fig. 4, the function R{a) decreases 
with increasing a. This means that when the particle is initially localized better than its Compton wavelength, the 
spreading of the packet occurs at a rate slower than in the non relativistic case. Moreover, it is easy to show that 
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in which the standard nonrelativistic result is corrected by the function F{a). Also this function decreases with 
increasing argument (see Fig. (4)), and the conclusion drawn for the spreading extends to the Hciscnbcrg principle. 

The role of the relativistic corrections in eqs. (3.8) and (3.10) can be understood as follows. The indetermination in 
position induces an indetermination in momentum and thus in velocity, which is in turn responsible for the relativistic 
corrections that we are dealing with. The contributions derived from these corrections are, however, of higher order 
in h, as easily inferred from the series expansion of the functions R{a) and F{a). For example, the function R{a) at 
lowest order in a reads 



R{a) 



1 



(3.11) 



This means that the ordinary nonrelativistic limit is obtained by just keeping the previous results at the lowest order 
in a2. 
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FIG. 4: Plot of the functions R{a) (blue) and F{a) (green) defined in the text. 

In the case of a particle subject to a linear potential 

H — c \J (? + — fx, 

the equations of motion are written as 



X = c 



P 



dt c2 ^ "2 

and, thus, for the position operator we obtain 



p,H 



x{t) = x(0) 
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\Jm?c^ + {tf + p(0))2 - \/m2c2 + p2(o) 



(3.12) 
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The evaluation of average values becomes more complicated, but still requires the use of fractional operators, we will 
not further dwell on this point but note that the method of the fractional derivatives offers a useful tool to treat 
problems involving even more complicated potentials. Also in this case we have assumed that the Schrodinger and 
Heisenberg pictures are linked by the usual unitary transformation (we have indeed evaluated the average value of a 
given operator as {X{t)) — exp(z Ht/h) X\'^)). Even though such an assumption sounds trivial, it is not taken for 
granted. In particular, when a potential is included, its mathematical legitimacy should be studied in a more rigorous 
way. 

IV. PARAMETRIZAZION OF A UNIT MATRIX AND FRACTIONAL OPERATORS 



In the previous section we have ignored the existence of another procedure to eliminate the square roots, which is 
essentially the Dirac parameterization method in terms of Clifford numbers [15]. The price to be paid is the necessity 
of introducing a matrix realization of these numbers, which inevitably increases the dimensionality of the problem. 

Before proceeding further, we it believe worthwhile to treat the Dirac method using a slightly unconventional 
procedure, according to which we will point out that this method is a tool to get the square root ofa2x2ora4x4 
unit matrix in terms of Pauli or Clifford numbers, respectively. 

The use of the Pauli numbers (matrices)'* 
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Prom now on, an underlined letter denotes a matrix and a product of matrix by a (pseudo)difTerential operator will be denoted by a hat 
over an underlined letter. 
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with 

{a„aj}=2Sijl [ai,a^] =2ieijk^k {i, j,k = 1,2,3) , (4.2) 

allows us to estabilish the following identity {v = {vi,V2,V3), Vi G C): 

N = a-v^ v^/l. (4.3) 
The anti-commutative properties of the Pauli numbers allows us also to prove that 

from which it is easy to show that 

e"^ = cosh(yiV)l + smh{yN)N. (4.5) 

As an example, let us apply the parameterization (4.3) to the case of the relativistic one-dimensional Schrodinger 
equation. In this case, the relevant Hamiltonian (see Eq. (2.1) written in terms of the usual coordinate x,t) can be 
written as 

H_ = aipc + a^mc^ (4.6) 

and the associated Schrodinger equation reads 

ihd,^=("l'' (4.7) 

where ^ is a two-component wave function. This equation should not be confused with the Pauli equation, but it 
should be understood as a two-dimensional Dirac equation with the components of the wave function referring to the 

positive and negative energy states. 

By using Eq. (4.5), the solution of Eq. (4.7) can be written as 

^{t) =exp(- ^Htj 1(0) 



cos 



that represents a kind of oscillation between the two states, with negative and positive energies. This aspect of the 
problem can be better understood from the Heisenberg equations of motion associated with the Hamiltonian (4.6), 
namely 

X = cai -T-P = ih-^ai = —2imc^a2- (4.9) 
at at at 

This equation should be completed with the Hamilton equation of motion for the a^. matrices that can be written in 
a vector form as follows: 

-^a = ^xa ^ = 2^(7r,0,l). (4.10) 

d t ~Ac 

This equation is reminiscent of the Bloch-type equation obtained by Feynman, Vernon, and Hellworth in their re- 
handling of the two-state Schrodinger equation [16]. In this case the two states are represented by the positive and 
negative energy solutions and the coupling between positive and negative states is realized by the matrix , whose 
time dependence is just an oscillating function between positive and negative values. In a forthcoming paper we will 
discuss more thoroughly the physical properties of eqs. (4.7) and (4.10). 

Analogous considerations can be developed for the pseudoheat equation (1.4), that, in a two-component form, reads 

9t^=~(,i (4.11) 
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where the negative and positive solutions are now interpreted as backward and forward heat fluxes. 
In the case of 4 x 4 matrices, the Dirac paranietrization applies: 



N = a■v + Pq = N^/l 
where a, (i are the Clifl[ord numbers (matrices) realized as 



{N = Vt^2 + g2 g C), 
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(j,fc = l,2,3) . 
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We stress that, if the only request is to write N as the product of a number by the square root of a unit matrix, 
other parametrizations, like 



a ■ V + a^u , 



(ai,a2 



or 



N = NI3, 



(4.15) 



(4.16) 



can also be used. The last parametrization has the advantage of being expressed in terms of a diagonal matrix. In 
some sense this is the essence of the so-called Foldy-Woythusen transformation [12], and will be discussed in the 
forthcoming section. 

In order to see whether this discussion has any relevance to the theory of fractional derivatives and/or pseudo- 
operators, we consider the following 4x4 operator: 

= {l-dl)l. (4.17) 

The square root of this operator can be written in terms of the Clifford matrices introduced before as follows: 

.1/2 



o 



ia- d + f3 
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5, (fc = 1,2,3), 



and therefore the pseudoheat equation (1.4) can be rewritten in the following form 

- 1/2 

dr^ = -O t 

where $ is a four-component function. By introducing the matrices 

7^=^afe (A; = 1,2, 3), 

and multiplying both sides of Eq. (4.19) by fi, one has (r' = — r) 



§_dr'^ = {if-d + l) 



(4.18) 

(4.19) 
(4.20) 

(4.21) 



that, as it will be discussed in the forthcoming section, is just a Dirac- like form. The last remark suggest that a link 
can be found between fractional calculus and Clifford algebras; further comments about this aspect will be presented 
in the next sections. 



V. DIRAC EQUATION 

In this section we will treat the Dirac equation and its implications for the fractional and pseudodifferential operators 
calculus. By assuming that the parametrization (4.12) can be applied also in then case in which v is replaced by a 
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differential operator, the Dirac equation for a free particle with mass m can be written in the Schrodinger form with 
a Hamiltonian given by 

= ca-p + (5m(? , (5.1) 
i.e., in terms of the normalized variables introduced in Eq. (2.1) 

idr"^ = (a,-w + Pj ^. (5.2) 

Moreover, by multiplying both sides of this equation for the matrix /3 it assumes the form 

ipdr^ = S^ (5.3) 

with 

S_ = f-n + 1 = - 7r2^, (5.4) 

where the analogy with the pseudo-heat equation (4.21) becomes even closer. The evolution operator associated with 
Eq. (5.2) is given by 

U_{t) = cos( ^1 + 7r2 t) 1 - i ^'"^ ^-^ ^ (g ■ f + /3~^) , (5.5) 

Vl + TT^ ^ 

and the corresponding solution 

*(?7,r) = Hr)^(?7) (5.6) 

shows all the relevant features of the evolution of a free Dirac particle, including the Zitterbewegung, if cj) contains 

positive and negative energy components. 

The Heisenberg equation for the position and momentum operators reads 

d - 



d7^ 



■i fj, a - TT + P = 5, OtTt = Q , (5.7) 

i.e., in this formalism the velocity operator is associated with the Clifford matrices. As for the matrices one has 

^a = 2i{^- aHj , (5.8) 
that allows to get the solutions for the position operator in the form 

^(r) = r}(0) + tttH-^ + '-HT^ (a(0) - % HT^^ {e-^"^ - l) . (5.9) 

The first two terms in this equation coincides with Eq. (3.3); the last is a further term accounting for the interference 
between particles with positive and negative energy (the Zitterbewegung term) and can be removed by their decoupling. 
Being the matrix /3 diagonal, the parametrization (4.16) automatically satisfies such a request, and yields 

Since the eigenvalues of /3 have the values ± 1, the positive and negative energy solutions have velocities in opposite 
directions. By integrating Eq. (5.10), since the momentum remains constant, we get 

^(t) = m + p J—^ T. (5.11) 

Vl + TT^ 

We close this section by addressing the problem of the existence of alternative forms for the Dirac equation depending 
on the chosen parametrization. In 1971, in a not widespread known paper [17], Dirac himself proposed a four- 
component equation, not containing negative energy solutions and valid for particles with integer spin. This equation, 
written in a form closely similar to the electron rclativistic case, contains a different realization of the matrices a^. 
and /3. Without entering into the details of the physical implications of this type of equation, which will be discussed 
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elsewhere, we consider the problem from the point of view that inspired this paper, by noting that for the product of 
the square root of a number the unit matrix also the following parametrization applies 



N = K- w + i6r 
where k-^ = —c^, ^2 = ^^i) % = ^^'^ 



{N = VtysTpT^ g c) . 



(5.12) 



/ 1 0\ 

1 
-10 

V -10 0/ 



These matrices satisfy the same anticommutation relation of Eq. (4.14), except that 6^ 
introduce the following equation 



(5.13) 



1. Accordingly we can 



(5.14) 



It resembles the ordinary Dirac equation, but the associated Hamiltonian is a non-Hermitian operator. Although it has 
been derived using the correct correspondence with the relativistic Schrodinger Hamiltonian, it cannot be exploited 
in unitary evolution processes and its physical properties are completely different from the ordinary Dirac equation 
(see [10] for further comments on the use of the matrices and 6 in the derivation of Dirac- like matrices). 
It is also interesting to note that 



is such that 



N = K- w + Sr 



N =sfu 



(5.15) 



(5.16) 



and therefore it can provide the square root of the matrix —1 or even of the null matrix (when w = r ). 



VI. CONCLUSIONS 



In this paper we have discussed a number of problems showing how different analytical methods, from apparently 
uncorrected fields, can be merged to provide interesting and useful tools. We have presented a general view on the 
theory of pseudo-operators and of fractional derivatives and on their implications for problems regarding the relativistic 
heat and the Dirac equations. The latter case has been discussed adopting a method involving the extraction of the 
square root of the unit matrix. We have touched many points and each of them would deserve a separate and deeper 
treatment. In particular, we have stressed that alternative parametrization may lead to a meaningful alternative. 
This observation is not entirely new and, apart from the already quoted Dirac paper [17], the same problem was 
addressed also in Ref. [18]. In this paper the authors insist on what they call the Dirac "dichotomy" and show that 
the alternative derivation, based on a 2 x 2 matrix representation, leads to a two-component equation of the type 
discussed in Sec. IV. 

In the paper we have not discussed the problem of the discrete symmetries of Eq. (4.7). For differential equations 
with fractional derivatives this problem deserves particular care [19]. We only remark that for parity (P) and time 
reversal (T), the corresponding operators can be constructed in terms of the Pauli matrices in analogy to the procedure 
followed in Ref . [20]. 

A further point, which has not been touched in the paper, is the nature of equations like the relativistic heat 
equation (1.4), generically treated as an evolution equation. The square root operator does not allow a classification 
in the usual sense (elliptic, parabolic, hyperbolic) (see Ref. [21]). We have avoided this point, since we have dealt 
with a genuine Cauchy problem, namely an initial value problem concerning the solution of a differential equation 
first order in time. The techniques associated with the evolution operator method are therefore sufficient, provided 
that the operator is well defined. The inclusion of boundary conditions would imply mathematical problems and 
technicalities beyond the scope of the paper. 

We close the paper briefly mentioning the application of this method to the case of the wave propagation. It is well 
known that the non-relativistic Schrodinger equation is used in classical optics to treat the paraxial wave propagation, 
namely the case in which the propagation occurs without significant variation of the field with respect to the axis 
of propagation [22] . The equation is obtained as a consequence of the application of the so-called Leontovich-Fock 
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approximation to the Helmholtz equation [23]. The Hamiltonian describing the propagation of an optical ray in a 
medium with an index of refraction n is^ 

H = - - p2 . (6.1) 

If the medium is homogeneous n can be considered constant; the variable p denotes the beam divergence and the 
use of the standard procedure allows us to replace it with the operator — i A Vg, where q is the transverse coordinate 
and A is the wavelength of the propagating radiation. The "quantization" procedure, involving A instead of h, is 
guaranteed by the fact that p and q are canonically conjugated variables. The method we have described in this paper 
can therefore be applied to the study of nonparaxial wave evolution and to the quantitative analysis of the deviation 
from the paraxial case. 



Appendix A 

The identity 



e 



is sometimes referred to as the Doetsch integral transform [7]. We dwell on the proof of this identity because, albeit 
available in literature [24, 25], it yields a first glimpse of how methods from different fields of calculus can be combined 
to get the desired result. 

The substitution f = 1/^^ in the integral in the right-hand side of Eq. (Al) leads us to consider the function defined 

as 

I{a, h) = ^ di exp ^- a" " |) > (A2) 
which is easily shown to satisfy the differential equation 

dbl{a,b) + 2 al{a,b) = Q, l{a,0) = ^ . (A3) 

The solution of this equation is given by 

I{a,b) = 7(a,0)e-2"^ , (A4) 
that for a = 1/2, b = x ^ reproduces Eq. (Al). In this case, from Eq. (A4) and (A2) we obtain also the identity 

that can be exploited to write the solution of Eq. (1.1) in the following alternative form 

e r2 



F{x,t) = ^ / d^ exp ( - ^ - ^ a, ) j{x) 



Jo 



^ ^'" d^e^^/Vfx- J) • (A6) 



^ This Hamiltonian is not invariant under Lorentz transformations and therefore cannot be considered relativistic. Nevertheless, a Dirac- 
like parametrization can be performed as 

-& = — a-p + i/3n. 



15 



Appendix B 

In this appendix we show how the aheady quoted Laplace transform method can be applied to evaluate the square 

root of a matrix. 

As an example, let us consider the following matrix 

R = al + ba^, (Bl) 

with a > 6 > 0. For its square root we can write 

^ Rf^rds'^), (B2) 



y/R Jo y/s 

where identity (2.4) has been used, and, as a consequence of Eq. (4.5), it turns out that 



-j= \ As — -j=— {cosh(&s)l — sinh(6 s)a;i} 



This method can be applied to more general matrix forms, by combining the previous integral transform with 
more conventional tools, like the Cayley-Hamilton theorem. Moreover, it is not limited to the square root but can be 
formulated in a wider context. For > 0, we get indeed 
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